rm(list=ls())

library(lme4)
library(readxl)
library(arm)
library(stargazer)

setwd("")

data_ab <- read_excel("data_abortion.xlsx")
data_ind <- read_excel("data_catindep.xlsx")

## ABORTION: 

# Factor DV and reference category, abortion
data_ab$language <- as.factor(data_ab$language)
data_ab$language <- relevel(data_ab$language, ref = "polish")
table(data_ab$cat_num_binary)

# Factor DV and reference category, Catalan independence
data_ind$language <- as.factor(data_ind$language)
data_ind$language <- relevel(data_ind$language, ref = "spanish")
table(data_ind$cat_num_binary, data_ind$language, data_ind$model)

## Model abortion
mlmod_ab <- glmer(cat_num_binary ~ language + (-1 + language|model), family=binomial("logit"), data=data_ab)
summary(mlmod_ab)
ranef(mlmod_ab)
fixef(mlmod_ab)
se.ranef(mlmod_ab)

## Model Catalan independence
mlmod_ind <- glmer(cat_num_binary ~ language + (-1 + language | model), family=binomial("logit"), data=data_ind)
summary(mlmod_ind)
ranef(mlmod_ind)
fixef(mlmod_ind)
se.ranef(mlmod_ind)

# Table 1 - Abortion (1)
row_english <- c(fixef(mlmod_ab)["languageenglish"], ranef(mlmod_ab)$model[,"languageenglish"])
row_english_se <- c(summary(mlmod_ab)["coefficients"][[1]]["languageenglish","Std. Error"], 
                    se.ranef(mlmod_ab)$model[,"languageenglish"])
row_swedish <- c(fixef(mlmod_ab)["languageswedish"], ranef(mlmod_ab)$model[,"languageswedish"])
row_swedish_se <- c(summary(mlmod_ab)["coefficients"][[1]]["languageswedish","Std. Error"], 
                    se.ranef(mlmod_ab)$model[,"languageswedish"])

# Table 1 - Catalan Idependence (2)
row_catalan <- c(fixef(mlmod_ind)["languagecatalan"], ranef(mlmod_ind)$model[,"languagecatalan"])
row_catalan_se <- c(summary(mlmod_ind)["coefficients"][[1]]["languagecatalan","Std. Error"], 
                    se.ranef(mlmod_ind)$model[,"languagecatalan"])

intercepts_row <- c(fixef(mlmod_ab)[1], fixef(mlmod_ind)[1])
se_intercepts_row <- c(summary(mlmod_ab)["coefficients"][[1]]["(Intercept)","Std. Error"],
                       summary(mlmod_ind)["coefficients"][[1]]["(Intercept)","Std. Error"])

# TABLE COMPONENTS:

# Coeffiecients and their standard errors:
rbind(row_english, row_english_se, row_swedish, row_swedish_se, row_catalan, row_catalan_se)
# Intercepts and their standard errors:
rbind(intercepts_row, se_intercepts_row)
# Number of observations in the data
c(nrow(data_ab), nrow(data_ind))
# Log likelihood
c(summary(mlmod_ab)["logLik"][[1]], summary(mlmod_ind)["logLik"][[1]])


## EXPORT TABLE: 

## Extract fixed & random effects for abortion model
fe_ab     <- summary(mlmod_ab)$coefficients
re_ab     <- ranef(mlmod_ab)$model
re_ab_se  <- se.ranef(mlmod_ab)$model

## Extract fixed & random effects for Catalan independence model
fe_ind     <- summary(mlmod_ind)$coefficients
re_ind     <- ranef(mlmod_ind)$model
re_ind_se  <- se.ranef(mlmod_ind)$model


add_rows <- function(est, se){
  r1 <- sprintf("%.3f", est)
  r2 <- sprintf("(%.3f)", se)
  return(rbind(r1, r2))
}

ab <- rbind(
  add_rows(fe_ab["languageenglish","Estimate"], fe_ab["languageenglish","Std. Error"]),
  add_rows(re_ab["GPT3","languageenglish"], re_ab_se["GPT3","languageenglish"]),
  add_rows(re_ab["GPT4","languageenglish"], re_ab_se["GPT4","languageenglish"]),
  
  add_rows(fe_ab["languageswedish","Estimate"], fe_ab["languageswedish","Std. Error"]),
  add_rows(re_ab["GPT3","languageswedish"], re_ab_se["GPT3","languageswedish"]),
  add_rows(re_ab["GPT4","languageswedish"], re_ab_se["GPT4","languageswedish"])
)
ind <- rbind(
  add_rows(fe_ind["languagecatalan","Estimate"], fe_ind["languagecatalan","Std. Error"]),
  add_rows(re_ind["GPT3","languagecatalan"], re_ind_se["GPT3","languagecatalan"]),
  add_rows(re_ind["GPT4","languagecatalan"], re_ind_se["GPT4","languagecatalan"])
)
int <- rbind(
  add_rows(fe_ab["(Intercept)","Estimate"], fe_ab["(Intercept)","Std. Error"]),
  add_rows(fe_ind["(Intercept)","Estimate"], fe_ind["(Intercept)","Std. Error"])
)

table <- rbind(
  English = c(ab[1,], ab[3,], ab[5,],"","",""),
  Blank1        = c(ab[2,], ab[4,], ab[6,],"","",""),
  Swedish = c(ab[7,], ab[9,], ab[11,],"","",""),
  Blank2        = c(ab[8,], ab[10], ab[12,],"","",""),
  Catalan = c("","","",ind[1,], ind[3], ind[5,]),
  Blank3      = c("","","",ind[2,], ind[4], ind[6,]),
  Intercept = c(int[1,], "", "", int[3,], "", ""),
  Blank4         = c(int[2,], "", "", int[4,], "", ""),
  Observations  = c("", 3000, "", "", 2000, ""),
  LogLik  = c("", -1678.912, "", "", -805.995, "")
)

rownames(table) <- c("English", " ", "Swedish", " ", "Catalan", " ", "Intercept", " ", "Observations", "Log Lik.")

stargazer(table,
          summary = FALSE,
          rownames = TRUE,
          out = "table1.html",
          type = "html")

stargazer(table,
          summary = FALSE,
          rownames = TRUE,
          out = "table1.tex",
          type = "latex")




